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TECHNICAL MEMORANDUM 1403 

ON THE INSTABILITY OF METHODS FOR THE INTEGRATION 
OF ORDINARY DIFFERENTIAL EQUATIONS 1 
By Heinz Rutishauser^ 


In spite of the remarkable publication of J. Todd (ref. 1) the 
essential points of which are related below, the author has since 
observed several times methods for the nummerical integration of differ- 
ential equations which, although subject only to a temptingly small 
truncation error, nevertheless involve the great danger of numerical 
instability. I want to state beforehand that this danger hardly exists 
for the well-tested methods of Runge-Kutta and Adams (extrapolation 
methods) if they are applied correctly. 

It is a natural characteristic that a differential equation to, be 
solved numerically is approximated by a difference equation, and that 
the latter is then solved. In order not to be_ forced to select an all 
too small interval, one prefers difference equations which approximate 
the differential equation as closely as possible but in compensation 
are of higher order than the original differential equation. Precisely 
in this, however, there lies a danger because the difference equation 
thereby has a greater diversity of possible solutions, and it may well 
happen that the numerical integration yields precisely one of the 
extraneous solutions which only at the beginning is in any way related 
to the desired solution of the differential equation. In the paper of 
J. Todd mentioned before several examples of this type have been 
enumerated. 

Consideration of the pertinent variation equation is particularly 
informative. It is very well possible that the differential variation 
equation is stable, that is, that it contains only converging solutions, 
whereas the difference variation equation is unstable since it possesses, 
due to the increased diversity of solutions, aside from the converging 
solutions, also solutions which increase exponentially. A deviation 
from the correct solution, once it exists, small as it may be - and such 


1 "Uber die Instabilitat von Methoden zur Integration gewShnlicher 
Dif f erentialgleichungen , " ZAMP, Kurze Mitteilungen, vol. Ill, 1952, 
pp. 65 - 74 . 

^Institut fur angewandte Mathematik der ETH, Zurich. 



2 


TM 1405 


deviations are unavoidable, because of the rounding-off errors - there- 
fore increases rapidly and may finally falsify considerably the solution 
obtained. Yet - we want to emphasize this once more - this instability 
is caused only by an inappropriate integration method. 

In the following discussion, several customary methods are examined 
from this viewpoint, and simple criteria for the stability of such meth- 
ods are indicated. For the rest, this report does not deal with error 
estimates. 


DIFFERENTIAL EQUATIONS OF THE FIRST ORDER 


y' = f (x,y) 


Variation equation 


dy 

(a) Integration by Means of Simpson T s Rule^ 

y k+i - y k-i + 1 ( y, k+i + 4y, k + y, k-i) (1 > 

This relationship, together with the differential equation, yields 
two equations for the unknown quantities y ^ +1 and y’^ which are 

solved mostly by iteration. If the differential equation is linear or 
quadratic in y, the iteration can be avoided. We assume, however, that 
one passes over, in any case, to the next integration interval only when 
the relationship (l) is satisfied. 

The difference variation equation pertaining to (l) is evidently 


^The phenomenon was observed on this example, and correctly inter- 
preted, also by Mr. G. Dahlquist, Stockholm (Lecture at the GaJMM- 
convention 1951 at Freiburg im Breisgau). Compare also: Z. angew. 

Math. Mech., 51 , 239, 1951. 

Sphere h signifies the length of the integration interval, and 
y^ stands as an abbreviation for y(kh) . 


* 
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^k+lP- “ j f y,k+l 



•y,k‘k 


" ( 1 + I f 7,k-]) Vl 


If one assumes f y to be constant, and chooses the expression T| k = 

for the solution of this equation, one obtains for X a quadratic 
equation with the solutions 


h 2 ^ 2 


hf,. 


= 1 + + "2" fy +•••»» e ^ 


"X — 1 4, h. f ll x* 2 

*2 “ - 1 + 3 f y " £g f y + * * * 




One recognizes easily that, of the two fundamental solutions ^ = ^1 
and T l2,k = ^2^ °f the difference variation equation, the first one 

approximates the solution of the differential variation equation whereas 
the second one is brought in by the numerical method. 


In particular 


~ (-dV^V 3 


represents for f y <0, thus precisely when the differential equation is 

stable, an oscillation which is slowly exponentially increasing. This 
has the effect that a small disturbance of the numerical solution - 
caused by a rounding-off or a truncation error - is intensified in the 
further course of the integration and finally gets completely out of 
hand. In Collatz 1 (ref. 2) book, the phenomenon is denoted as "roughening 
phenomenon"; means for elimination of this inconvenience are given. On 
the other hand, the explanation given there is not complete; the phe- 
nomenon occurs only for fy < 0; there is nothing to be apprehensive of 

for f y ^ 0 which is very important particularly in regard to the ordi- 
nary Simpson’s integration rule ^f y = 0 j . 


(b) Integration According to Runge-Kutta and Similar Methods 


Since these methods calculate y^+l ^ rom y^. according to a 
prescribed rule and without use of the preceding values y^-^, y^_g, . . 
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the order remains unchanged in the transition from the differential 
equation to the difference equation; thus no foreign solutions are hrought 
in, and no instability is to he feared. 

The same property can be found in a method indicated by W. E. Milne 
(ref. 3). 


(c) Integration According to Adams 
We consider a four -point formula 

^k+i - + (sy'k+i + - sr'k-i + y‘k-2^ + h5 • • • (2) 

This again yields, together with the differential equation, two equa- 
tions for the unknown quantities y k+1 and y*k+i are mostly 

solved by iteration. 

The difference variation equation pertaining to ( 2 ) becomes 

^ f y,k+l^k+l - ^ ^ f y,k^k + §£ ^y, k-l^k-1 “ f y,k-2 t, k-2 = 0 

If one again considers fy as constant, the expression = A^ yields 
an equation of the third order for A. A solution An of this equation 

hfy k X 

lies very close to e , therefore k = ^1 corresponds to the 

* k 

solution of the differential variation equation whereas 7^ and 7^ 
are extraneous solutions. 

However, the equation for A is reduced toA^ - A^ = 0 when h 
tends toward 0 so that, for a sufficiently small h, one will have at 

any rate small A2 and Aj, namely A/ ±\J-hfy/ 24 . The extraneous solu- 

k Tr 

tions T] 2 ,k = A2 and 'Hj^k = ^3 thus converge rapidly. For a suffi- 
ciently small h, Adams' method is therefore stable. 

(d) Variants of (c) 

In order to improve the accuracy of Adams' method, one may use also * 

other expressions for the corrector instead of ( 2 ) . As long as the 
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corresponding five- or six-point formulas are involved, there are no 
objections, but one has to be careful when y^+l n °4 calculated 

from y k and the derivatives as in (2) but perhaps from y k _ 2 . or 

y k _5 and the derivatives, as for instance in 

yk+1 = yk-3 + H^y'k+1 + 32y'k + ^'k-l + 32y‘ k _2 + Ty'k-3^ + h? * * * 

In fact the pertinent difference variation equation has a solution 
T\ k = with A = - Jl - £ y + . . 

so that the method is unstable for fy < 0. 

DIFFERENTIAL EQUATIONS OF THE SECOND ORDER 
y" = f(x,y,y’) 

Insofar as these equations are solved by separation into a system 
of two equations of the first order, what was said so far is valid. 
Particularly in the case of numerical integration of damped oscillations 
we must caution against the methods (a) and (d) . 

However, there exist also methods which solve an equation of the 
second order without transformation into a system: 


(e) The Method of Central Differences? 

The formulas on which this method is based are (especially for 
second order) 

y k +i = 2yk - ^k-i + + 10 y"k + y"k-^ 


(4) 


y‘k + i - y, k-i + |(y"k + i + 4 y"k + y Vi) 


(5) 


^Compare reference 2, p. 80. 
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They yield, together with the differential equation, three equations 
for the unknown quantities y k+1 , y ' k+1 , and y" k+1 * The two simul- 
taneous difference variation equations pertinent to (4) and ( 5 ) are 
solved with the. expression T) fc = pX , T)’^ = qX , under the assumption 

of a constant f v and f 1 ; because of 

J J 


I’m - f yVl + 


one obtains, with the abbreviations a for h2fy/l2 and b for 
hfyi/5, the equations ' 

p[h 2 - 2X + 1) - a(X 2 + 10X + l)J - 
q | b(X 2 + 10X + 1) = 0 [from (4)] 

q£(X 2 - 1) - b(X 2 + 4X + 1)J - p ^ a(X 2 + 4X + l) =0 j^from (5)J 


( 6 ) 


These equations can exist simultaneously with (p,q) (0,0) only when 

the determinant of this equation system for p and q vanishes; one 
obtains after some calculations 

(X - l) j^(X^ - l) (X - l) - a(X + 1) (X^ + 10X + l) - 
b(X - l)(X^ + 4X + 1)J = 0 

The four solutions of these equations are 

where cc.q and Og are the solutions of 


Xq = 1 + aj h + . . 
Xg — 1 + agh + . . 

x 5 = 1 

^ - 5 f y' + • 


the equation a?- - ctfy t - fy = 0 
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Evidently A^ and A^ are the regular solutions of the difference 

variation equation; they correspond to two fundamental solutions of 
the differential variation equation; A^ and in contrast, are 

extraneous. As long as ^y« = 0 , there is nothing to fear, in partic- 
ular, the method may he strongly recommended for a y'-free equation, 
hut for fy» < 0 (damped oscillations) 

%,k “ \ k ~ (-l) k e~( kh /3)f y t 

increases, and Aj, too, may still become dangerous because A^ ^ = 1 

also becomes finally very large, compared to a function converging to- 
ward zero. 

The author completely calculated the example y" + y* + 1 . 25 y = 0 
with the initial conditions y = 0 , y* = 1 (exact solution: 

f ) 0 

e sin x/ on the sequence- controlled computing machine of the ETH. 

There follow a few excerpts from the thus obtained table of func- 
tions (we calculated with h = 0 . 1 ) : 


X 

y 

y' 

4.8 

-0.0903699 

0.0531227 

4.9 

- .0847792 

.0584842 

5.0 

- .0787132 

.0626410 

5-1 

- .0722891 

.0656573 

5.2 

- .0656173 

.0676070 


In this region nothing conspicuous is noticeable yet, the y-values 
deviate from the true values approximately by one in the last decimal 
place, and only formation of the differences for the y' -values reveals 
a certain irregularity. For t = 17 , however, the influence of A^ k 

becomes pronounced for the y* -values and also for the differences of 
the y-values: 
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X 

y 

y' 

17.0 

17.1 

17.2 

17.3 

17.4 

17.5 

-0.00019574 

- .00019061 

- .00018366 

- .00017524 

- .00016548 

- .00015475 

0.00005017 

.00005253 

.00008620 

.00008239 

.00011235 

.00010258 


The considerably weaker oscillation of the y- values follows also 
from the equations (6): for A = A^ one obtains from the first of 

these equations 

= - — f, r i , here therefore p ~ 

q. 4 4 6 y 6o 0 

The further course of numerical integration does not require any comment;: 


X 

y 

y’ 

22.8 

-0.00000815 

0.00005320 

22.9 

- .00000864 

- .00006078 

23.0 

- .00000862 

.00005968 

23.1 

- .00000887 

- .00006247 

23.2 

- .00000861 

.00006601 

23.3 

- .00000868 

-.00006486 

29.5 

- .00000140 

-.00053037 

29.6 

.00000041 

.00054889 

29.7 

- .00000144 

- .OOO56693 

29.8 

.00000050 

.00058682 

29.9 

- .00000148 

.00060603 

30.0 

.00000060 

- .00062735 


The author is well aware that the assumption of a ponstant fy and 
f v i in the above considerations greatly restricts the generality. How- 

t7 

ever, the results show that what matters is only the sign of these quan- 
tities, and this sign is indeed invariable in a great many cases. The 
statements are, therefore, qualitatively almost generally valid. Only 
when fy changes its sign, from time to time in the course of the inte- 
gration, for instance when method (a) is being used, a special case 
arises since the occurring oscillations alternately increase and are 
damped again. 
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GENERAL CONSIDERATIONS 




■» 


Consideration, of the examples suggests the conjecture that insta- 
bility may occur precisely in the case of integration methods which 
form. y k+1 by integration of y* over several intervals (two intervals 

in method (a), one interval in (b) and (c), four intervals in (d) , two 
intervals in (e)). However, this is not exactly the case and we shall 
therefore subject a general integration method to an examination. 

Almost all known methods use relationships which are contained in 
the general formula 


N 



. -1 
k+1 



(n-1) (n-l) 
n-l,rk+j 



a ( n_1 ) 

nj 



+ 


+ 


. N-n+1 
h 


■^=- (n-1) 

Z ®Nj 

-m 
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There are, in addition, differential equations 


F n ^k+l' y k+l' * • *> y k+l) = 0 . 

F n+i ( x k+i' y k+i' • • v ar£j 1} ) * 0 ^ 


( 8 ) 




M X k+l' y k+l’ 



) 


which form, together with the equations (7), N + 1 equations for the 

(N) 

N + 1 unknowns y k+1 , y*^^, . . ., y^ + ^. Generally, N = n; however, 

W. E. Milne '(ref. 3 ) uses in the method previously mentioned higher 
derivatives than those appearing in the differential equation. There- 
fore he differentiates them several times in order to obtain the required 
number of relationships. One may thus obtain the equations F n+1 * * * 

to Fjy by differentiating the initial equation F n . 

If one, furthermore, brings everything in the equations (j) to one 
side, the variation equations for the entire system (j) and (8) read 

evidently (with a^ = -lj 


N 


SI SI 

(iai j=-m 


(i). u-i (u) 
a v /h p ti vh / = 0 
k+j 


(i “ 0, 1, . • ., n — i) 



dFj, 

dy(n) 



= 0 


(i a n, n + 1, . . ., N) 


( Lll * 

If one uses for this the expression tk = pA J , one obtains in sub- 

J 

stitutlng a system of N + 1 simultaneous equations for the p^ which 
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can "be satisfied only when the determinant of the system disappears 


N / 1 


5 £ - 0 


(1 — 0, lj • • • y n — 1} 


dF l 

£l >> P,1 = 0 


(i = n, n + 1, . . N) 


,(i) 


If one defines, in addition, with the coefficients a^.' appearing in 
the formulas (7) the functions 


VW ■ t ^ 


J—m 


follows 


D(7v) = 


III 

=i 

<r 

0 for | 

AOO 

0 

• 

• 

hAoi 

All 

• 

9 

• 

• 

• 

0 

0 

• • • • 

F n,y 

v • • • 

F n,y‘ 

F n+l,y 

• 

F n+l,y' 

• 

• 

%,y 

F N,y' 


Ai-l^n-l 


b? "lA] w 


t-N-n+l* , mj 

n A n-1,N 


r n)Jr (i) o o 


‘ P N, y (M) 



12 


TM 1403 


If the method Is to reproduce exactly a polynomial of the ith degree, 
which is the solution of a differential equation, together with its 
derivatives - and this one may require - the conditions 


y(i) = 1 y(i+l) S y(i+2) = y(i+3) s ... = y(N) = 0 
J 3 A 0 3 


must be compatible. Hence follows, however, by substitution into the 
ith of the equations (7) : V a^ 1 ) = 0, therefore A ii (l) = 0. 

J J 


The equation D(A) =0 which is decisive for the stability of the 
method must have n solutions in the neighborhood of X = 1, corre- 
sponding to the n independent solutions of the variation equation. 

In fact one finds for h = 0 where D(A) is, except for one factor, 
reduced to AqcAu * * * -^n-l n-1 ' fcha ' fc A = 1 is an n-fold zero of 

D(A), because of A^Cl) = 0. 


All other zeros of D(A) correspond therefore to extraneous solu- 
tions of the difference equation; in order to make the method stable, 
they must lie, for a sufficiently small h, in the interior or at most 
on the periphery of the unit- circle. This is certainly the case when 
for h a 0 all zeros of D(A) lie in the interior of the unit circle, 
and certainly not when individual ones are outside it. Therefore: 

Sufficient condition for the stability of the method ( 7 ) for 
sufficiently small h: 


All functions 


Aii(A) 



a 


(i)^j+m 


possess, aside from the trivial simple zero A = 1 only zeros with 

| a l < 1. 

Necessary condition : None of the functions Aj^(A) has a zero 

Outside the unit circle. 

. Todd considered methods which satisfy not even the necessary 
condition. The solution obtained then becomes completely useless 
already after a few intervals. 
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The remaining functions 


when the necessary condition, 
satisfied. 


A i[ ,(A) can influence the stability only 
but not the sufficient condition, is 


APPLICATIONS 


For the formula (l) there results (one has N = n * 1, m =» 1) 


Aqo = -A 2 + 1 


Aq! - 1 (A 2 + 4A + 1) 
5 


F, «= y* - f(x,y) 


Therefore 


D(?0 = 


1 - A c 


-tv 


(A 2 + k-K + 1) 
1 


The fact that A qq has two zeros with |a| a 1 already suggests cau- 
tion, but moreover one reads off imediately that L(A) is positive 
for f y < 0 and A = -1, and negative. In contrast, for A ■ -«>. 

Thus a 2ero lies to the left of -1; the method Is unstable. 

For the formula 5*^-2 in the book of Collatz mentioned (p. 81), 
there is (N = n =» k, m = l) 


= y 17 - f(x,y,y , ,y",y m ) 
AqQ = App = -(A - l) 2 


A- 


'•ll 

12 


Ajj = -A + 1 

2A 

= i (A 2 + + 1) 


All other A^ occur only with at least h 2 in the determinant, 
one has, except for terms with h 2 


Thus 
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D(A) = 


■ (A - l) 2 
0 
0 
0 

-f\r 


0 

-A 2 + 1 
0 
0 


-f. 


y* 


0 

2Ah 

■(A - l) 2 
0 


-v 


0 

0 

0 

-A 2 + 1 


-f 


y"‘ 


0 
0 
0 

k (A 2 + 4A + 1)1 


Far A = D is positive, for A = -1 - e, D has the sign of 


-A 2 + 1 | (A 2 + 4A + 1) 

5 

-v 1 

■Which for fy.ni < 0 and a sufficiently small e is negative . Therefore 

this method is unstable for f m < 0. 

«/ 


On the other hand it is easy to indicate methods which are always 
stable. One need only shape the formulas ( 7 ) in such a manner that 
every line begins with 


•y 


M . y (i) 


k+1 




a (i) v (i+l) 

a i+l,j yV j+k ; 


+ h 


(i = 0,1, 


n-l) 


Thereby A^(A) = -A IIH '^ + A m and has therefore only the trivial zero 
A = 1 on the periphery of the unit circle. 


SUMMARY 


In the numerical solution of a differential equation as a difference 
equation, the latter is usually of higher order and therefore has more 
solutions than the original differential equation. It may well be that 
some of these "extra" solutions grow faster than any solution of the 
given equation; in this case the computational solution has the tend- 
ency to follow one of these and has after a certain number of integra- 
tion steps nothing to do with the original differential equation. 
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The author gives some examples and a criterion for stability of 
integration methods. This criterion is then applied to some well-known 
integration formulas. 


Translated by Mary L. Mahler 
National Advisory Committee 
for Aeronautics 
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